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ABSTRACT 

A viscoelastic internal variable constitutive theory is applied to a higher-order elastic 
beam theory and finite element formulation. The behavior of the viscous material in the 
beam is approximately modeled as a Maxwell solid. The finite element formulation 
requires additional sets of nodal variables for each relaxation time constant needed by the 
Maxwell solid. Recent developments in modeling viscoelastic material behavior with strain 
variables that are conjugate to the elastic strain measures are combined with advances in 
modeling through-the-thickness stresses and strains in thick beams. The result is a viscous 
thick-beam finite element that possesses superior characteristics for transient analysis since 
its nodal viscous forces are not linearly dependent on the nodal velocities, which is the case 
when damping matrices are used. Instead, the nodal viscous forces are directly dependent 
on the material's relaxation spectrum and the history of the nodal variables through a 
differential form of the constitutive law for a Maxwell solid. The thick beam quasistatic 
analysis is explored herein as a first step towards developing more complex viscoelastic 
models for thick plates and shells, and for dynamic analyses. 

The internal variable constitutive theory is derived directly from the Boltzmann 
superposition theorem. The mechanical strains and the conjugate internal strains are shown 
to be related through a system of first-order, ordinary differential equations. The total time- 
dependent stress is the superposition of its elastic and viscous components. Equations of 
motion for the solid are derived from the virtual work principle using the total time- 
dependent stress. Numerical examples for the problems of relaxation, creep, and cyclic 
creep are carried out for a beam made from an orthotropic Maxwell solid. 

INTRODUCTION 

Advanced composites technology has led to the use of highly viscous, low-modulus 
materials that are combined with the traditional high-modulus, load carrying materials to 
produce stiff, highly damped structures. The quasi-static and dynamic analyses of such 
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structures require improvements in the material damping representation over the standard 
proportional damping schemes. Halpin and Pagano 1 demonstrated that the relaxation 
moduli for anisotropic solids produce symmetric matrices that can be expanded in a Prony 
series form (i.e., a series of exponentially decaying terms). Early viscoelastic models for 
small deformations of composites focused on computing the complex moduli for 
anisotropic solids from the elastic properties of the fibers and the complex modulus 
properties of the matrix 2,3 . Recently, various classical constitutive models have been used 
including generalized Maxwell and Kelvin-Voigt solids 4,5 . These constitutive models have 
practical value since they provide adequate approximations for the dynamic softening and 
hysteresis effects - the phenomena that are not directly proportional to strain rates. 

In this paper, a brief review of the history integral form of the Maxwell solid is 
presented to provide background for the differential constitutive law. The interested reader 
is referred to Coleman and Noll 6 and Schapery 7 for more comprehensive discussions on the 
subject. A new differential constitutive law for the Maxwell solid is then derived. This 
form is a special case of the model developed by Johnson and Stacer 8 for large strain 
viscoelastic deformations of rubber. It has also been used by Johnson et al. 9 ' 10 to formulate 
a viscoelastic, large-displacement shell finite element. The differential constitutive law is 
then combined with the higher-order beam theory and finite element formulation of 
Tessler 11 providing viscoelastic capability for thick beams. The additional strain variables 
required in the constitutive law are replaced with element nodal variables that are conjugate 
to the elastic nodal variables. This results in only minor modifications to the elastic finite 
element code. Finally, several numerical examples are presented demonstrating the 
viscoelastic, thick-beam response under quasi-static loading. 

MAXWELL SOLID IN HISTORY INTEGRAL FORM 

The stress-strain relation for a linear elastic material can be written in the tensor form as 

°ij ~ Cijkl £ kl (1) 

where G are the stress components, C are the elastic stiffness coefficients, and £ kl 
are the strains. When a linear viscoelastic material is subjected to an instantaneous 
incremental strain, A£ w , the time dependent stress takes the form 

< J ij(t) = C ijkl A£ kl 


( 2 ) 
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where the viscous stresses, oJj(t), are monotonically decreasing functions of time. The 
Boltzmann superposition method is often used to approximate Eq(2) as follows. The 
viscous stresses, <T - (/), are factored such that 

G Jj (0 ~ Cijki ( 3 ) 

The functions (t) are referred to as time dependent moduli. These monotonically 
decreasing functions are approximated in time using a Prony series, i.e., 

N — — 

C\ ijkl (0 = Cjjkl X e " ( 4 ) 

n=l 

where T„, C* w > 0. The stresses then become 

N — — 

G ij (0 = Cijki^ £ ki + X Cijki e T " Ae kl (5) 

n= 1 

The above approximation is extended to the case of a continuously deforming solid by 
associating the continuous time dependent strain with an incremental strain history and 
convoluting Eq(5) in time. The approximation to the time dependent stresses becomes 

M N M _( , Z bn) 

°ij (0 = Cijki 'L A m £ ki + X C *jki X H(t -t m )e T " A m £ kl (6) 

m = 1 n = 1 m = 1 

where the strain increments are set at times t m for m = 1 , M , and use is made of the 
Heaviside unit step function, H(t — t m ) . At this juncture, it is often customary to define 
the viscous moduli as functions of the relative time, t — t m , i.e., 

N 

Cijki ( ? - ) = X Cijki ~ t m) e n C) 

n = 1 

The constitutive model in Eq(6) then takes the form 

M M 

G ij (0 = C ijkl ^ A m e u + ^ Cjj kl ( t - t m )A m e kl (8) 

m - 1 m — 1 

Assuming that strains are smooth functions of time, and taking the limit as 
A m t = t m+ j — t m — > 0 for all m, gives rise to 
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where it is noted once again that the viscous moduli, Cijtdi 1 ~ T )’ are monotonically 

decreasing in time. Eq(9) is often referred to as the history integral for a linear Maxwell 
solid. 

In many practical applications, adequate time-dependent stress predictions can be 
obtained with only several terms in the Prony series, Eq(4). The constitutive model as 

presented in Eq(9) requires that the history of the measurable kinematics, £ kl (T), be 
known in addition to the Prony series. This leads to computational algorithms that must 
determine how much of the history to retain in order to update the viscous stress 
approximation as time evolves. 

MAXWELL SOLID IN DIFFERENTIAL FORM 

Following Johnson et al., 8 10 new constitutive equations for a linear viscous solid are 
derived. The new constitutive equations are in differential form and they are equivalent to 
the history integral form just described. Departing from the history integral formulation at 

Eq(6), defining internal strain variables, *£ w , which relate to the strains as, 

A mn£kl(t-*m) = H(t-t m )e ^ A m £ kl (10) 

introducing Eq(10) into Eq(6) and factoring out C* Jkl , the stresses appear as 

M N M 

® ij (0 — Cykl ^ ^m^kl Cjjkl ^ ^ A m n £ k i {t ~ t m ) ( 11 ) 

m = 1 n=l m - 1 

Following the Maxwell solid formulation, it is assumed that the strains are smooth 
functions of time. In the limit as A m t = t m+l —t m — > 0 for all m, Eq(l 1) becomes 

N t 

<*ij (0 = C m e u {t) + C* jkl X J d*£ kl (t- r) ( 12 ) 

n = 1 x=-°° 

It is desirable to derive a differential equation for the time dependent strain variables. In the 
limit as A m t 0, Eq(10) becomes 

(f-T) 

= tn d£ kl {T) 

Integrating Eq(13) with respect to the history, T, yields 


( 13 ) 
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„*£«(')= J 


r ^ k ,(r) 


Xn dr 


Differentiating Eq(14) with respect to the current time, t , yields 


d* n £ki( 0 


J_ r jgafr) 

r J /'/'r 


r " dr + 


(f) 


Substituting Eq(14) into Eq(15) yields the differential equations for the internal strain 
variables in the form 


d n £ kl 


for each n . 


Introducing Eq(14) into Eq(12) results in the stress equation given by 


a ij (0 - Cijkl £ kl ( f ) + Cjjkl X , n £ kl (0 


Eqs(16) and (17) represent the constitutive model in differential form. Note that this 
particular form is for the case of a material with a relaxation modulus given by Eq(4). 
Also, for a material whose modulus is expressed by Eq(4), the constitutive model of 
Eqs(16) and (17) is equivalent to the history integral model given by Eq(9). In what 
follows, the differential form of the constitutive model is explored in the context of a 
higher-order beam theory and its associated finite element. 

VISCOELASTIC HIGHER-ORDER BEAM 

In this section a quasi-static Maxwell solid version of Tessler’s" higher-order beam 
theory is formulated and a simple beam finite element is derived. The beam dimensions and 
sign convention are shown in Figure 1 . The viscoelastic constitutive model for the beam 
that is consistent with Eqs(16) and (17) can be written in matrix form as 


s M = Ce+*cjr *e(() 


7 * * 


for each n. 


where 
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S r =(< 7 «. °ZZ' T xz) 

® — ’ ^zz * y xz ) 

* T ( * * * 

n® *“(n^xc» n^zz* nY xz 



Ql 

Q 3 

0 " 


~*C„ 

* c 

'-'13 

0 
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C13 

C33 

0 
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*c 13 
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'-'33 
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0 

0 

Cs5_ 


0 

0 

v 

'-'55 


In this higher-order theory, the components of the displacement vector are approximated 
through the beam thickness by way of five kinematic variables, i.e., 

u x (x,z,t) = u(x,t) + hC,0(x,t) (20) 

U z (x,z,t) = w(x,t) + %w x (xj) + ^ 2 -ijw 2 (*,/) (21) 

where £ = z! h denotes a nondimensional thickness coordinate and 2 h is the total 


thickness. Note that u(x, t) represents the midplane (i.e. reference plane) axial 
displacement, 9(x, t ) is the bending rotation of the cross-section of the beam, Vf(jc, t) is 
the weighted-average deflection," and W, (x, t) and w 2 (x,t ) are the higher-order 
transverse displacement variables enabling a parabolic distribution of H, (jc, Z, t) through 
the thickness. The above displacement assumptions give rise to the following thickness 
distributions for the strains: a linear axial strain, a cubic transverse normal strain, and a 
quadratic transverse shear strain. 1 1 These strain components have the following form 


= u(x,t), x +h£6(x,t), x 

(22) 


(23) 

= <t>xz( 0 +0(x,O) 

(24) 


where 


<h(0= 


hv l3 C{4-7i 2 ) 
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. U0 = 


14fA(3-f 2 ) 
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* 0 = 




The simplest finite element approximation of this beam theory, as explored in Ref. [11], 
involves a three-node configuration (see Figure 2) which is achieved by the following 
interpolations 


U (T], t) = (1 - ri)u 0 ( t ) + T\u x (f) 


(25) 
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9(ri, t) = (1 — 77)0 O (f) + 770 ] (7) (26) 

0 = (! - 0 V 0 (0 + 0 W 1 (0 - ^ 0(! - l){Oo (0 - #i (0) ( 27 ) 

W] (ri,t) = W ] (t) (28) 

w 2 {n^) = W 2 (0 (29) 


where Tf = x / 1 is the nondimensional element length coordinate. Note that the nodal 
degrees-of-freedom at the two ends of the element are subscripted with indices 0 and 1 . 
Since the strains do not possess derivatives of the w^(T],t) and W 2 (77, t) variables, 
these variables need not be continuous at the element nodes and, hence, their simplest 
approximation is constant for each element. Their corresponding degrees-of-freedom are 
attributed to a node at the element midspan. 

For a quasi-static loading, the virtual work statement for an element of volume V 
with the differential form of the Maxwell constitutive law included can be written as 


JVc&k/V + XJ *e T *C8* n edV-8W = 0 


(30) 


n = 1 


where the first integral represents the internal virtual work done by the elastic stresses, the 

second is the internal virtual work done by the viscous stresses, and 8W is the virtual 
work done by the external forces. Introducing Eqs(25)-(29) into Eqs(20)-(21) and 
substituting the results into Eqs(22)-(24) yields finite element approximations of the strains 
in terms of the nodal variables, i.e., 


e = Bu 


(31) 


where 


B = 


1 

1 

0 

0 - 


0 

0 

K 
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_ z 
i 

i 
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j_ 

h 
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1 

1 

^ 0 


0 

0 

h 2 
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0 - 
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z 

1 

*x_ 

1 

txz 


(32) 


and U T = [uq , vv 0 , 0 q , Wj , W 2 , u x , W] , 0] ) denotes the element nodal displacement 
vector. 


^ $ 

A set of analogous nodal variables, n U, and corresponding viscous strains, „ e , are 
introduced. These are related by 
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nC = B W U (33) 

The virtual work statement for an element then becomes 

N 

u 7 ’jB r CBJVAi+X« u7 J Br * CB ^n u -^ =0 (34) 

n = 1 

By defining the integrals in Eq(34) as stiffness matrices , there results 
N 

u r k &i+^ *u 7 *k <5*u- 8W =0 ( 35 ) 

n = 1 

Since Eq(19) implies <5e —S n t when t is constant, the virtual work takes on a simpler 
form 


u r k 


N 


+s> 

n = 1 


T * 1 


da -8W =0 


This implies that at any given time the element equilibrium equations are 

N 

k U + y k*u = f for each element 

n = 1 


(36) 


(37) 


where f denotes the element consistent load vector due to the external loading. 
Introducing Eqs(31) and (33) into the differential equations for the strain variables in 
Eq(19) yields 


d*U *11 da 

1 — — for each n (38) 

dt Z n dt 

The global equilibrium equations are determined by the standard assembly of the 


element equations, Eqs(37). Note, there is no assembly for Eqs(38). The variables * U 
are independent from element to element (recall, these variables carry the time dependent 
information for the material within the element). The global equilibrium equations at a 
given time are 

g — ^ mech — ^vi'.vr (39) 


where 11 g denotes the global nodal variable vector, K is the elastic stiffness matrix, 
F mech * s global force vector due to mechanical loads, and is the assembled 

vector for ^ k n U . The viscoelastic problem is solved by simultaneously integrating 
n = 1 
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the differential and algebraic equations expressed by Eqs(38) and (39), where the latter is 
subject to the appropriate boundary restraints. 

As far as the finite element implementation is concerned, a conventional linearly elastic 
code can be readily adapted to perform the viscoelastic analysis for a Maxwell material, 
i.e., for a material whose relaxation stiffness coefficients can be modeled with a Prony 

series. First, the instantaneous stiffness coefficients, *C J y W , are used in place of the elastic 

values to compute the element viscous stiffness matrices, k, which are stored for 

* 

repeated use. The internal nodal variables for each element, n U, are set equal to their 
initial values and stored. A predictor-corrector algorithm is then used to integrate Eqs(38) 
and (39) in time. The predictor-corrector integration algorithm used in this effort is 
described in the Appendix. 


APPLICATIONS 

Numerical solutions representative of stress relaxation, creep, and cyclic creep for a 
thick orthotropic beam are presented. The beam elastic stiffness coefficients ( C matrix in 
Eq (18)) for the state of plane stress can be written in terms of engineering material 
constants as 


Qi - 


1 — V V 

xz zx 


C 33 - 


1 — v v 

Y xz zx 


* C n —v xz C 33 , C , 55 -G 


55 


xz 


A unidirectional E-glass/epoxy laminate is considered for which the material constants are: 

E x =3S.6GPa, E z = 8.27 GPa , G xz =4.14 GPa, v xz =0.26, and 

V zx = V XZ E X / E, . The viscous relaxation properties were computed from complex 
modulus vs. frequency data for the E-glass/epoxy reported in Ref[12]. The equations for 
the real and imaginary parts of the modulus of a Maxwell series were least-squares fit to the 
data in a frequency range of 45 Hz - 145 Hz. The least squares fit was performed with the 
constraint that the moduli in the Maxwell series each be positive. The series was defined 
with ten time constants; T n = 1.0E-4, 3.162E-4, 1.0E-3, 3.162E-3, 1.0E-2, ... , 1.0E+1, 
and infinity. The Prony series was scaled so that its equivalent static value (at t equal to 
infinity) was unity. The resulting series is 

P(t) = 1.0 + 0.01755 tT 00001 ' + 0.000257 e~° mt + 0.072014 <r° 3162 ' 


The time dependent stiffness values for the E-glass/epoxy are given by C v = C P(t ) . 

The beam dimensions are as follows: L = 0.1m, 2h = 0.02m, andb = 0.01m (refer 
to Figure 1). 



10 


Example 1 . A cantilever beam with w, u, 9 fixed at point A has a prescribed 
deflection w at point B that is ramped from 0 to -1 cm in 0.05 sec and then held constant 
at -1 cm. Figure 3 depicts the value of the maximum axial stress computed at point D as a 
function of time. Also shown are the elastic and viscous stress components comprising the 
total stress. The decay of the total viscoelastic stress to its elastic value as time is increased 
demonstrates the expected step-strain relaxation behavior. 

Example 2 . A simply-supported beam with w and U fixed at point A and w fixed at 

point B is subject to a uniform, top-surface pressure, q + (f). The time-dependent value of 
the pressure is ramped from 0.0 to 1 .0 MPa in 0.05 sec and then held constant. Figure 4 
shows the value of the W deflection at the midspan of the beam. The viscous solution 
shows the expected creep response. 

Example 3 . The simply-supported beam in the preceding example is subject to a 
harmonic pressure loading given by 

q + (t) = 0, t < 0 and q + (t) = 0.05 * (l - cos(507Ef)) MPa, t> 0 

Figure 5 depicts the transverse normal stress and strain values at the midspan (point C) 
versus time. The upward drift of the maximum and minimum values of the strain 
demonstrates cyclic creep. 

CONCLUSION 

A differential form of the Maxwell viscous solid constitutive theory has been derived 
and implemented within a higher-order-theory beam finite element. The finite element 
formulation is attractive for several reasons: (1) The constitutive constants are the same as 
those needed in the classical history-integral model, and they are also readily available from 
step-strain relaxation tests, (2) The state variables are conjugate to the elastic strain 
measures; hence, they are consistent with the kinematic assumptions of the elastic 
formulation, (3) The update of the state variables can be performed in a parallel computing 
environment, allowing the viscous force vector in the equations of motion to be determined 
efficiently within the predictor-corrector algorithm, (4) Applications of time-dependent 
displacements and loads are performed within the same finite element algorithm, and (5) 
The higher-order beam theory accounts for both transverse shear and transverse normal 
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deformations — the effects that need to be accounted for in thick and highly orthotropic 
beams and high-frequency dynamics. 

The computational examples for the problems of relaxation, creep, and cyclic creep 
clearly demonstrated the predictive capabilities of the finite element formulation. 
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Stresses and displacements computed at C & D. 

Figure 1. Thick-beam geometry, kinematics, and loading. 



Figure 2. A three-node, higher-order-theory thick-beam element. 
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Figure 3. Cantilevered beam under prescribed tip deflection. 
Stress in top surface at support. 



Figure 4. Simply-supported beam under prescribed pressure loading. 
Weighted-average deflection, w, at center of beam. 
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Figure 5. Simply-supported beam under harmonic pressure loading. 

Transverse normal strain and stress at beam center versus time. 


APPENDIX 

The predictor-corrector integration algorithm is described below. The storage 
requirements, beyond the requirements for the elastic problem, involves two sets of 

vectors, * U, for each element, two global vectors, U g ,a global force vector, F y(SC , a full 

set of element viscous stiffness matrices, *k, and a matrix which carries the Prony series 
information for each element. The time integration algorithm used in this effort provides a 
trapezoidal solution to the differential equations while simultaneously solving the algebraic 
equilibrium equations. The elastic (static) finite element code already contains the assembly 
and solver subroutines needed. Building the viscous (quasi-static) finite element code only 
requires that additional storage be allotted and a few modified call statements be added. 
The algorithm is outlined below. 

Time Integration Algorithm 

A. Initialize variables { \ =$ t ,2 => t + At) 
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* u 2 = 0 V elements. 
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S> 2 


= 0 


P\n = 


( x _ J*_\ 

2t„ 


1 + 


Ar 


2r 


integration constant. 


« J 


Pi n ~ 


1 + 


At 


It 


integration constant. 


n ) 


B. Move data to next time step. 
*il[ = *ll 2 Velements. 

p = F 

r mech y \ r mech, 2 

F r: F 

r vise,] r vise ,2 


U S,1 = U £.2 

t = / + A t 

^ mech,2 ~ ^ mech (0 

U g,2 = K ‘ [ ¥ mech, 2 ~ F wc,l ] seed for Ste P C 

C. Update internal variables (element components of global 
vectors implied). 

>2 =A„>1 +P2n( u g,2 

D. Update viscous forces. 

N 

¥ vise, 2 assemble £ *k*U 2 
n = 1 

E. Update global displacements. 

U S,2 =K~ l [ ¥ mech,2 ~ ¥ visc,2 ] 

F. If changes to 2 are small then go to Step B. 

Else go to Step C. 
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